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Abstract 

Hydraulic fracturing of deep shale beds to develop natural gas has caused concern regarding the potential for 
various forms of water pollution. Two potential pathways — advective transport through bulk media and preferential 
flow through fractures — could allow the transport of contaminants from the fractured shale to aquifers. There 
is substantial geologic evidence that natural vertical flow drives contaminants, mostly brine, to near the surface 
from deep evaporite sources. Interpretative modeling shows that advective transport could require up to tens of 
thousands of years to move contaminants to the surface, but also that tracking the shale could reduce that transport 
time to tens or hundreds of years. Conductive faults or fracture zones, as found throughout the Marcellus shale 
region, could reduce the travel time further. Injection of up to 15,000,000 L of fluid into the shale generates 
high pressure at the well, which decreases with distance from the well and with time after injection as the fluid 
advects through the shale. The advection displaces native fluids, mostly brine, and fractures the bulk media 
widening existing fractures. Simulated pressure returns to pre-injection levels in about 300 d. The overall system 
requires from 3 to 6 years to reach a new equilibrium reflecting the significant changes caused by tracking the 
shale, which could allow advective transport to aquifers in less than 10 years. The rapid expansion of hydraulic 
fracturing requires that monitoring systems be employed to track the movement of contaminants and that gas 
wells have a reasonable offset from faults. 



Introduction 

The use of natural gas (NG) in the United States has 
been increasing, with 53% of new electricity generating 
capacity between 2007 and 2030 projected to be with NG- 
fired plants (EIA 2009). Unconventional sources account 
for a significant proportion of the new NG available to 
the plants. A specific unconventional source has been 
deep shale-bed NG, including the Marcellus shale primar- 
ily in New York, Pennsylvania, Ohio, and West Virginia 
(Soeder 2010), which has seen over 4000 wells devel- 
oped between 2009 and 2010 in Pennsylvania (Figure 1). 
Unconventional shale-bed NG differs from conventional 



Hydrologic Consultant, 6320 Walnut Creek Road, Reno, 
NV 89523; (775) 530-1483; fax: (775) 530-1483; tom_ 
myers@charter.net 

Received August 201 1 , accepted February 201 2. 

© 2012, The Author(s) 

Ground Water © 201 2, National Ground Water Association, 
doi: 10.1 1 1 1/j'.1 745-6584.201 2.00933.x 



sources in that the host-formation permeability is so low 
that gas does not naturally flow in timeframes suitable for 
development. Hydraulic fracturing (fracking, the industry 
term for the operation; Kramer 2011) loosens the forma- 
tion to release the gas and provide pathways for it to move 
to a well. 

Fracking injects up to 17 million liters of fluid 
consisting of water and additives, including benzene at 
concentrations up to 560 ppm (Jehn 2011), at pressures 
up to 69,000 kPa (PADEP 2011) into low permeability 
shale to force open and connect the fractures. This is 
often done using horizontal drilling through the middle 
of the shale with wells more than a kilometer long. The 
amount of injected fluid that returns to the ground surface 
after fracking ranges from 9% to 34% of the injected fluid 
(Alleman 2011; NYDEC 2009), although some would be 
formation water. 

Many agency reports and legal citations (DiGiulio 
etal. 2011; PADEP 2009; ODNR 2008) and peer- 
reviewed articles (Osborn et al. 2011; White and Mathes 
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Figure 1. Location of Marcellus shale in the northeastern United States. Location of Marcellus wells (dots) drilled from July 
2009 to June 2010 and total Marcellus shale wells in New York and West Virginia. There are 4064 wells shown in Pennsylvania, 
48 wells in New York, and 1421 wells in West Virginia. Faulting in the area is documented by PBTGS (2001), Isachsen and 
McKendree (1977), and WVGES (2011, 2010a, 2010b). 



2006) have found more gas in water wells near areas 
being developed for unconventional NG, documenting the 
source can be difficult. One reason for the difficulty is 
the different sources; thermogenic gas is formed by com- 
pression and heat at depth and bacteriogenic gas is formed 
by bacteria breaking down organic material (Schoell 
1980). The source can be distinguished based on both 
C and H isotopes and the ratio of methane to higher chain 
gases (Osborn and Mcintosh 2010; Breen et al. 2007). 
Thermogenic gas can reach aquifers only by leaking from 
the well bore or by seeping vertically from the source. 
In either case, the gas must flow through potentially very 
thick sequences of sedimentary rock to reach the aquifers. 
Many studies which have found thermogenic gas in water 
wells found more gas near fracture zones (DiGiulio et al. 
2011; Osborn et al. 2011; Breen et al. 2007), suggesting 
that fractures are pathways for gas transport. 

A pathway for gas would also be a pathway for flu- 
ids and contaminants to advect from the fractured shale to 
the surface, although the transport time would be longer. 
Fracking fluid has been found in aquifers (DiGiulio et al. 



2011; EPA 1987), although the exact source and pathways 
had not been determined. With the increasing development 
of unconventional NG sources, the risk to aquifers could 
be increasing. With so little data concerning the movement 
of contaminants along pathways from depth, either from 
wellbores or from deep formations, to aquifers, conceptual 
analyses are an alternative means to consider the risks. 

The intent of this study is to characterize the risk 
factors associated with vertical contaminant transport 
from the shale to near-surface aquifers through natural 
pathways. I consider first the potential pathways for 
contaminant transport through bedrock and the necessary 
conditions for such transport to occur. Second, I estimate 
contaminant travel times through the potential pathways, 
with a bound on these estimates based on formation 
hydrologic parameters, using interpretative MODFLOW- 
2000 (Harbaugh et al. 2000) computations. The modeling 
does not, and cannot, account for all of the complexities 
of the geology, which could either increase or decrease 
the travel times compared to those considered herein. 
The article also does not include improperly abandoned 
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boreholes which could cause rapid transport in addition 
to natural pathways. 



Method of Analysis 

Using the Marcellus shale region of southern New 
York (Figure 1), I consider several potential scenarios 
of transport from shale, 1500 m below ground surface 
(bgs) to the surface, beginning with pre-development 
steady state conditions to establish a baseline and 
then scenarios considering transport after fracking has 
potentially caused contaminants to reach formations 
above the shale. To develop the conceptual models and 
MODFLOW-2000 simulations, it is necessary first to 
consider the hydrogeology of the shale and the details 
of hydraulic fracturing, including details of how fracking 
changes the shale hydrogeologic properties. 

Hydrogeology of Marcellus Shale 

Shale is a mudstone, a sedimentary rock consisting 
primarily of clay- and silt- sized particles. It forms 
through the deposition of fine particles in a low energy 
environment, such as a lake- or seabed. The Marcellus 
shale formed in very deep offshore conditions during 
Devonian time (Harper 1999) where only the finest 
particles had remained suspended. The depth to the 
Marcellus shale varies to as much as 3000 m in parts 
of Pennsylvania, and averages about 1500 m in southern 
New York (Soeder 2010). Between the shale and the 
ground surface are layers of sedimentary rock, including 
sandstone, siltstone, and shale (NYDEC 2009). 

Marcellus shale has very low natural intrinsic perme- 
ability, on the order of 10 -16 Darcies (Kwon et al. 2004a, 
2004b; Neuzil 1986, 1994). Schulze-Makuch et al. (1999) 
described Devonian shale of the Appalachian Basin, of 
which the Marcellus is a major part, as containing "coaly 
organic material and appear either gray or black" and 
being "composed mainly of tiny quartz grains < 0.005 mm 
diameter with sheets of thin clay flakes." Median particle 
size is 0.0069 ± 0.00141 mm with a grain size distribu- 
tion of <2% sand, 73% silt, and 25% clay. Primary pores 
are typically 5 x 10 -5 mm in diameter, matrix porosity 
is typically 1% to 4.5% and fracture porosity is typically 
7.8% to 9% (Schulze-Makuch et al. 1999 and references 
therein). 

Porous flow in unfractured shale is negligible due 
to the low bulk media permeability, but at larger scales 
fractures control and may allow significant flow. The Mar- 
cellus shale is fractured by faulting and contains synclines 
and anticlines that cause tension cracks (Engelder et al. 
2009; Nickelsen 1986). It is sufficiently fractured in some 
places to support water wells just 6 to 10 km from where 
it is being developed for NG at 2000 m bgs (Loyd and 
Cars well 1981). Conductivity scale dependency (Schulze- 
Makuch et al. 1999) may be described as follows: 

K = Cv m 



K is hydraulic conductivity (m/s), C is the intercept of a 
log-log plot of observed K to scale (the K at a sample 
volume of 1 m 3 ), v is sample volume (m 3 ), and m is 
a scaling exponent determined with log-log regression; 
for Devonian shale, C equals 10 -14 3 , representing the 
intercept, and m equals 1.08 (Schulze-Makuch et al. 
1999). The very low intercept value is a statistical but 
not geologic outlier because it corresponds with very 
low permeability values and demonstrates the importance 
of fracture flow in the system (Schulze-Makuch et al. 
1999). Most of their 89 samples were small because the 
deep shale is not easily tested at a field-scale and no 
groundwater models have been calibrated for flow through 
the Marcellus shale. Considering a 1-km square area with 
30-m thickness, the Kh would equal 5.96 x 10 -7 m/s 
(0.0515 m/d). This effective K is low and the shale would 
be an aquitard, but a leaky one. 

Contaminant Pathways from Shale to the Surface 

Thermogenic NG found in near- surface water wells 
(Osborn et al. 2011; Breen et al. 2007) demonstrates the 
potential for vertical transport of gas from depth. Osborn 
et al. (2011) found systematic circumstantial evidence for 
higher methane concentrations in wells within 1 km of 
Marcellus shale gas wells. Potential pathways include 
advective transport through sedimentary rock, fractures 
and faults, and abandoned wells or open boreholes. Gas 
movement through fractures depends on fracture width 
(Etiope and Martinelli 2002) and is a primary concern for 
many projects, including carbon sequestration (Annunzi- 
atellis et al. 2008) and NG storage (Breen et al. 2007). 
Open boreholes and improperly sealed water and gas 
wells can be highly conductive pathways among aquifers 
(Lacombe et al. 1995; Silliman and Higgins 1990). 

Pathways for gas suggest pathways for fluids and 
contaminants, if there is a gradient. Vertical hydraulic 
gradients of a up to a few percent, or about 30 m over 
1500 m, exist throughout the Marcellus shale region as 
may be seen in various geothermal developments in 
New York (TAL 1981). Brine more than a thousand 
meters above their evaporite source (Dresel and Rose 
2010) is evidence of upward movement from depth to 
the surface. The Marcellus shale, with salinity as high 
as 350,000 mg/L (Soeder 2010; NYDEC 2009), may 
be a primary brine source. Relatively uniform brine 
concentrations over large areas (Williams et al. 1998) 
suggest widespread advective transport. The transition 
from brine to freshwater suggests a long-term equilibrium 
between the upward movement of brine and downward 
movement of freshwater. Faults, which occur throughout 
the Marcellus shale region (Figure 1) (Gold 1999), could 
provide pathways (Konikow 2011; Caine et al. 1996) 
for more concentrated advective and dispersive transport. 
Brine concentrating in faults or anticline zones reflects 
potential preferential pathways (Wunsch 2011; Dresel and 
Rose 2010; Williams 2010; Williams et al. 1998). 

In addition to the natural gradient, buoyancy would 
provide an additional initial upward push. At TDS equal 
to 350,000 mg/L, the density at 25 °C is approximately 
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1290 kg/m 3 , or more than 29% higher than freshwater. 
The upward force would equal the difference in weight 
between the injected fluid and displaced brine. As an 
example, if 10,000,000 L does not return to the surface as 
flowback (Jehn 2011), the difference in mass between the 
volume of fracking fluid and displaced brine is approxi- 
mately 3,000,000 kg, which would cause an initial upward 
force. The density difference would dissipate as the salt 
concentration in the fracking fluid increases due to diffu- 
sion across the boundary between the fluid and the brine. 

In just Pennsylvania, more than 180,000 wells had 
been drilled prior to any requirement for documenting 
their location (Davies 2011), therefore the location of 
many wells is unknown and some have probably been 
improperly abandoned. These pathways connect aquifers 
through otherwise continuous aquitards; overpressuriza- 
tion of lower aquifers due to injection near the well 
pathway could cause rapid transport to higher aquifers 
(Lacombe et al. 1995). In the short fracking period, the 
region that is overpressurized remains relatively close to 
the gas well (see modeling analysis below), therefore it 
should be possible for the driller to locate nearby aban- 
doned wells that could be affected by fracking. This article 
does not consider the potential contamination although 
unlocated abandoned wells of all types must be considered 
a potential and possibly faster source for contamination 
due to fracking. 

Effect of Hydraulic Fracturing on Shale 

Fracking increases the permeability of the targeted 
shale to make extraction of NG economically efficient 
(Engelder et al. 2009; Arthur et al. 2008). Fracking 
creates fracture pathways with up to 9.2 million square 
meters of surface area in the shale accessible to a 
horizontal well (King 2010; King et al. 2008) and 
connects natural fractures (Engelder et al. 2009; King 
et al. 2008). No post-fracking studies that documented 
hydrologic properties were found while researching this 
article (there is a lack of information about pre- and post- 
fracking properties; Schweitzer and Bilgesu 2009), but 
it is reasonable to assume the K increases significantly 
because of the newly created and widened fractures. 

Fully developed shale typically has wells spaced at 
about 300-m intervals (Edwards and Weisset 201 1; Soeder 
2010). Up to eight wells may be drilled from a single 
well pad (NYDEC 2009; Arthur et al. 2008), although 
not in a perfect spoke pattern. Reducing by half the 
effective spacing did not enhance overall productivity 
(Edwards and Weisset 2011) which indicates that 300-m 
spacing creates sufficient overlap among fractured zones 
to assure adequate gas drainage. The properties controlling 
groundwater flow would therefore be affected over a large 
area, not just at a single horizontal well or set of wells 
emanating from a single well pad. 

Fracking is not intended to affect surrounding forma- 
tions, but shale properties vary over short ranges (King 
2010; Boyer et al. 2006) and out-of-formation fracking is 
not uncommon. In the Marcellus shale, out-of-formation 
fracks have been documented 500 m above the top of the 



shale (Fisher and Warpinski 2011). These fractures could 
contact higher conductivity sandstone, natural fractures, or 
unplugged abandoned wells above the target shale. Also, 
fluids could reach surrounding formations just because of 
the volume injected into the shale, which must displace 
natural fluid, such as the existing brine in the shale. 

Analysis of Potential Transport along Pathways 

Fracking could cause contaminants to reach overlying 
formations either by fracking out of formation, connecting 
fractures in the shale to overlying bedrock, or by 
simple displacement of fluids from the shale into the 
overburden. Advective transport, considered as simple 
particle velocity, will manifest if there is a significant 
vertical component to the regional hydraulic gradient. 

Numerical modeling, completed with the MODF 
LOW-2000 code (Harbaugh et al. 2000), provides flex- 
ibility to consider potential conceptual flow scenarios, but 
should be considered interpretative (Hill and Tiedeman 
2007). The simulation considers the rate of vertical trans- 
port of contaminants to near the surface for the different 
conceptual models, based on an expected, simplified, real- 
istic range of hydrogeologic aquifer parameters. 

MODFLOW-2000 is a versatile numerical modeling 
code, but there is insufficient data regarding the geology 
and water chemistry between aquifers and the deep shale, 
such as salinity profiles or data concerning mixing of the 
brine with fracking fluid, to best use its capabilities. As 
more data becomes available, it may be useful to consider 
simulating the added upward force caused by the brine by 
using the SEAWAT-2000 module (Langevin et al. 2003). 

Vertical flow would be perpendicular to the general 
tendency for sedimentary layers to have higher horizontal 
than vertical conductivity. Fractures and improperly 
abandoned wells would provide pathways for much 
quicker vertical transport than general advective transport. 
This article considers the fractures as vertical columns 
with model cells having much higher conductivity than 
the surrounding bedrock. The cell discretization is fine, so 
the simulated width of the fracture zones is realistic. Dual 
porosity modeling (Shoemaker et al. 2008) is not justified 
because turbulent vertical flow through the fractures is 
unlikely, except possibly during the actual fracking that 
causes out-of-formation fractures, a scenario not simulated 
here. MODFLOW-2000 has a module, MNW (Halford 
and Hanson 2002), that could simulate rapid transport 
through open bore holes. MNW should be used in 
situations where open boreholes or improperly abandoned 
wells are known or postulated to exist. 

The thickness of the formations and fault would affect 
the simulation, but much less than the several-order-of- 
magnitude variation possible in the shale properties. The 
overburden and shale thickness were set equal to 1500 and 
30 m, respectively, similar to that observed in southern 
New York. The estimated travel times are proportional 
for thicker or thinner sections. The overburden could 
be predominantly sandstone, with sections of shale, 
mudstone, and limestone. The vertical fault is assumed 
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to be 6-m thick. The fault is an attempt at considering 
fracture flow, but the simulation treats the 6-m wide fault 
zone as homogeneous, which could underestimate the real 
transport rate in fracture-controlled systems which could 
be highly affected by dispersion. The simulation also 
ignores diffusion between the fracture and the adjacent 
shale matrix (Konikow 2011). 

There are five conceptual models of flow and trans- 
port of natural and post-fracking transport from the level 
of the Marcellus shale to the near- surface to consider 
herein: 

1. The natural upward advective flow due to a head 
drop of 30 m from below the Marcellus shale to the 
ground surface, considering the variability in both shale 
and overburden K. This is a steady state solution for 
upward advection through a 30-m thick shale zone 
and 1500-m overburden. Table 1 shows the chosen K 
values for shale and sandstone. 

2. Same as number 1, but with a vertical fracture 
connecting the shale with the surface, created using 
a high-conductivity zone in a row of cells extending 
through all from above the shale to the surface. This 
emulates the conceptual model postulated for flow into 
the alluvial aquifers near stream channels, the location 
of which may be controlled by faults (Williams et al. 
1998). The fault K varies from 10 to 1000 times the 
surrounding bulk sandstone K (K ss ). 

3. This scenario tests the effect of extensive fracturing 
in the Marcellus shale by increasing the shale K 
(K S h) from 10 to 1000 times its native value over 
an extensive area. This transient solution starts with 
initial conditions being a steady state solution from 
scenario 1. The K s ^ increases from 10 to 1000 times 
at the beginning of the simulation, to represent the 
relatively instantaneous change on the regional shale 
hydrogeology imposed by the fracking. The simulation 
estimates both the changes in flux and the time for the 
system to reach equilibrium. 

4. As number 3, considering the effect of the same 
changes in shale properties but with a fault as in 
number 2. 

5. This scenario simulates the actual injection of 13 to 
17 million liters of fluid in 5 d into fractured shale 
from a horizontal well with and without a fault. 

Model Setup 

The model domain was 150 rows and columns spaced 
at 3 m to form a 450-m square (Figure 2) with 50 layers 
bounded with no flow boundaries. The 30-m thick shale 
was divided into 10 equal thickness layers from layer 40 
to 49. The overburden layer thickness varied from 3 m 
just above the shale to layer 34, 6 m from layer 33 to 29, 
9 m from layer 28 to 26, 18 m in layer 25, 30 m from 
layer 24 to 17, 60 m from layer 16 to 6, 90 m from layer 
5 to 3, and 100 m in layers 2 and 1. A 6-m wide column 
from layer 39 to the surface is added for some scenarios 
in the center two rows to simulate a higher K fault. 



Table 1 

Sandstone (ss) and Shale (sh) Conductivity (K) 
(m/d) and the Steady State Flux (m 3 /d) for Model 
1 Scenarios 



Flux K ss K sh 



1.7 


0.1 


0.00001 


1.8 


0.5 


0.00001 


1.9 


1 


0.00001 


1.9 


5 


0.00001 


2.0 


10 


0.00001 


2.0 


50 


0.00001 


2.0 


100 


0.00001 


1.7 


0.1 


0.00001 


9.5 


0.1 


0.00005 


19.0 


0.1 


0.0001 


81.2 


0.1 


0.0005 


135.9 


0.1 


0.001 


291.5 


0.1 


0.005 


340.9 


0.1 


0.01 


394.3 


0.1 


0.05 


401.8 


0.1 


0.1 


409.2 


0.1 


0.5 


40.7 


0.001 


0.1 


186.0 


0.005 


0.1 


339.1 


0.01 


0.1 


988.3 


0.05 


0.1 


1297.3 


0.1 


0.1 


1748.0 


0.5 


0.1 


1826.1 


1 


0.1 


1902.8 


5 


0.1 


1915.4 


10 


0.1 


338.3 


0.1 


0.01 


984.1 


0.5 


0.01 


1292.5 


1 


0.01 


1731.5 


5 


0.01 


1816.0 


10 


0.01 


17.4 




0.0001 


86.3 




0.0005 


176.7 




0.001 


775.1 




0.005 


1292.5 




0.01 


2746.8 




0.05 


3183.2 




0.1 


3650.5 




0.5 


3719.9 




1 



The model simulated vertical flow between constant 
head boundaries in layers 50 and 1, as a source and 
sink, so that the overburden and shale properties control 
the flow. The head in layers 50 and 1 was 1580 and 
1550 m, respectively, to create a gradient of 0.019 over 
the profile. Varying the gradient would have much less 
effect on transport than changing K over several orders 
of magnitude and was therefore not done. 

Scenario 5 simulates injection using a WELL bound- 
ary in layer 44, essentially the middle of the shale, from 
columns 25 to 125 (Figure 2). It injects 15 million liters 
over one 5-d stress period, or 3030 m 3 /d into 101 model 
cells at the WELL. The modeled was changed to its 
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Figure 3. Sensitivity of the modeled head response to the 
storage coefficient used in the fractured shale for model layer 
39 just above the shale. 



Figure 2. Model grid through layer 44 showing the horizon- 
tal injection WELL (red) and DRAIN cells (yellow) used to 
simulate flowback. There is only one monitoring well because 
the off-center well is not used in layer 44. 



assumed fracked value at the beginning of the simula- 
tion. Simulating high rate injection generates very high 
heads in the model domain, similar to that found sim- 
ulating oil discharging from the well in the Deepwater 
Horizon crisis (Hsieh 2011) and water quality changes 
caused by underground coal gasification (Contractor and 
El-Didy 1989). DRAIN boundaries on both sides of the 
WELL simulated return flow for 60 d after the completion 
of (Figure 2), after which the DRAIN was deactivated. 
The 60 d were broken into four stress periods, 1, 3, 6, and 
50 d long, to simulate the changing heads and flow rates. 
DRAIN conductance was calibrated so that 20% of the 
injected volume returned within 60 d to emulate standard 
industry practice (Alleman 2011; NYDEC 2009). Recov- 
ery, continuing relaxation of the head at the well and the 
adjustment of the head distribution around the domain, 
occurred during the sixth period which lasted for 36,500 d. 

There is no literature guidance to a preferred value 
for fractured shale storage coefficient, so I estimated S 
with a sensitivity analysis using scenario 3. With fractured 
K S h equal to 0.001 m/d, two orders of magnitude higher 
than the in situ value, the time to equilibrium resulting 
from simulation tests of three fractured shale storage 
coefficients, 10 -3 , 10 -5 , and 10 _7 /m, varied twofold 
(Figure 3). The slowest time to equilibrium was for S = 
10 _3 /m (Figure 3), which was chosen for the transient 
simulations because more water would be stored in the 
shale and flow above the shale would change the least. 

Results 

Scenario 1 

Table 1 shows the conductivity and flux values 
for various scenarios. The steady state travel time 



for a particle through 1500 m of sandstone and shale 
equilibrates with one of the formations controlling the 
advection (Figure 4). For example, when the K s ^ equals 
1 x 10 -5 m/d, transport time does not vary with K ss . For 
K ss at 0.1 m/d, transport time for varying K s ^ ranges from 
40,000 to 160 years. The lower travel time estimate is for 
similar to that found by Schulze-Makuch et al. (1999). 
The shortest simulated transport time of about 20 years 
results from both the sandstone and shale K equaling 
1 m/d. Other sensitivity scenarios emphasize the control 
exhibited by one of the media (Figure 4). If is low, 
travel time is very long and not sensitive to K ss . 

Scenario 2 

The addition of a fault with K one to two orders of 
magnitude more conductive than the surrounding sand- 
stone increased the particle travel rate by about 10 times 
(compare Figure 5 with Figure 4). The fault K controlled 
the transport rate for ^ S h less than 0.01 m/d. A highly 
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Figure 4. Sensitivity of particle transport time over 
1500 m for varying shale and sandstone vertical K. 
Effective porosity equals 0.1. (1) — varying K ss , K S h = 
10" 5 m/d; (2)— varying K ss , K sh = 0.1 m/d; (3)— varying 
K ss , K sh = 0.1 m/d; (4)— varying K ss , K sh = 0.01 m/d; and 
(5)— varying K ss , K sh = 1.0 m/d. 



6 T.Myers GROUNDWATER 



NGWA.org 



« 3Ll .(rl 



10.QO0 ■ 



■con 



* 100 





0.COD001 QQQ031 Q.D001 



■ (1) 



o.ooi o.oi oi 

Conductivity {m/d) 



i ml; 



Figure 5. Variability of transport through various scenarios 
of changing the K for the fault or shale. Effective porosity 
equals 0.1. (1) — varying /f S h, # S h = 0.01 m/d; (2) — varying 
K S h, K S h = 0.1 m/d; (3) — no fault; (4) — varying K fault, 
K S h = 0.1 m/d, K S h = 0.01 m/d. Unless specified, the vertical 
fault has K = 1 m/d for variable K s ^. 



conductive fault could transport fluids to the surface in 
as little as a year for equal to 0.01 m/d (Figure 5). 
However, a fault did not significantly change the overall 
model flux, so with fault values are not shown in Table 1 . 

Scenarios 3 and 4 

Scenarios 3 and 4 estimate the time to establish a 
new equilibrium once the K s ^ changes, due to fracking, 
between values specified in scenarios 1 and 2. Equilibrium 
times vary by model layer as the changes propagate 
through the domain, and flux rate for the simulated 
changes imposed on natural background conditions. The 
fracking-induced changes cause a significant decrease in 
the head drop across the shale and the time for adjustment 
of the potentiometric surface to a new steady state depends 
on the new shale properties. 

The time to equilibrium for one scenario 3 simulation, 

changing from 10 -5 to 10 -2 m/d with K ss equal 
to 0.1 m/d, varied from 5.5 to 6.5 years, depending 
on model layer (Figure 6). Near the shale (layers 39 
and 40), the potentiometric surface increased from 23 
to 25 m reflecting the decreased head drop across the 
shale. One hundred meters higher, in layer 20, the 
potentiometric surface increased about 20 m. Simulation 
of scenario 4, with a fault with K = 1 m/d, decreased 
the time to equilibrium to from 3 to 6 years within the 
fault zone, depending on model layer (Figure 6). Highly 
fractured sandstone would allow more vertical transport, 
but advective flow would also increase so that the base 
K ss would control the overall rate. 

The flux across the upper boundary changed within 
100 years for scenario 3 from 1.7 to 345 m 3 /d, or 
0.000008 to 0.0017 m/d, reflecting control by K ss . There 
is little difference in the equilibrium fluxes between 
scenario 3 and 4 indicating that the fault primarily affects 
the time to equilibrium rather than the long-term flow rate. 
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Figure 6. Monitoring well water levels for specified model 
layers due to fracking of the shale; monitor well in the center 
of the domain, including in the fault, K of the shale changes 
from 0.00001 to 0.01 m/d at the beginning of the simulation. 



Scenario 5: Simulation of Injection 

The injection scenarios simulate 15 million liters 
entering the domain at the horizontal well and the 
subsequent potentiometric surface and flux changes 
throughout. The highest potentiometric surface increases 
(highest injection pressure) occurred at the end of injec- 
tion (Figure 7), with a 2400 m increase at the horizontal 
well. The simulated peak pressure both decreased and 
occurred longer after the cessation of injection with dis- 
tance from the well (Figures 7 and 8). The pressure at 
the well returned to within 4 m of pre-injection levels in 
about 300 d (Figure 7). After injection ceases, the peak 
pressure simulated further from the well occurs longer 
from the time of cessation, which indicates there is a pres- 
sure divide beyond which fluid continues to flow away 
from the well bore while within which the fluid flows 
toward the well bore. The simulated head returned to 
near pre-injection levels slower with distance from the 
well (Figure 7), with levels at the edge of the shale (layer 
40) and in the near-shale sandstone (layer 39) requir- 
ing several hundred days to recover. After recovering 
from injection, the potentiometric surface above the shale 
increased in response to flux through the shale adjusting to 
the change in shale properties (Figure 8), as simulated in 
scenario 3. The scenario required about 6000 d (16 years) 
for the potentiometric surface to stabilize at new, higher, 
levels (Figure 8). Removing the fault from the simulation 
had little effect on the time to stabilization, and is not 
shown. 

Prior to injection, the steady flux for in situ shale 
= 10 -5 m/d) was generally less than 2 m 3 /d and 
varied little with K ss (Figure 4). Once the shale was 
fractured, the sandstone controlled the flux which ranges 
from 38 to 135 m 3 /d as K ss ranges from 0.01 to 0.1 m/d 
(Figure 9), resulting in particle travel times of 2390 and 
616 years, respectively. More conductive shale would 
allow faster transport (Figure 4). Adding a fault to the 
scenario with K ss equal to 0.01 m/d increased the flux to 
approximately 63 m 3 /d and decreased the particle travel 
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Figure 7. Simulated potentiometric surface changes by layer 
for specified injection and media properties. The monitoring 
point is in the center of the domain. Fault is included. lf S h = 
0.01 m/d, K sh = 0.001 m/d. 5 (fractured shale) = 0.001/m, 
S (ss) = 0.0001/m. 
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Figure 9. Comparison of flux for three scenarios. Flowback 
is the same for all scenarios. (1): Kss = 0.01 m/d, Ksh = 
0.001 m/d, Fault K = 1 m/d; (2): Kss = 0.01 m/d, Ksh = 
0.001 m/d, no fault; (3) Kss = 0.1 m/d, Ksh = 0.001 m/d, no 
fault. 
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Figure 8. Simulated potentiometric surface changes for lay- 
ers within the shale and sandstone. CW is center moni- 
toring well and EW is east monitoring well, about 120 m 
from the centerline. Fault is included. The line for layer 
2, CW plots beneath the line for layer 2, EW. K ss = 
0.01 m/d, K sh = 0.001 m/d, S (fractured shale) = 0.001/m, 
S (ss) = 0.0001/m. 



time to 31 years. Approximately, 36 m 3 /d flowed through 
the fault (Figure 9). The fault properties control the 
particle travel time, especially if the fault K is two or 
more orders of magnitude higher than the sandstone. 

Simulated flowback varied little with ^ S h because it 
had been calibrated to be 20% of the injection volume. 
A lower storage coefficient or higher K would allow the 
injected fluid to move further from the well, which would 
lead to less flowback. 

Vertical flux through the overall section with a fault 
varies significantly with time, due to the adjustments in 
potentiometric surface. One day after injection, vertical 
flux exceeds significantly the pre-injection flux about 
200 m above the shale (Figure 10). After 600 d, the 
vertical flux near the shale is about 68 m 3 /d and in 



ico 




200 400 600 800 1000 1 
M«t#ft above Shale 
■ One toy after Injection — - 600 D»yi - - - 100 Yea^ 



Figure 10. Upward flux across the domain section as a 
function of distance above the top of the shale layer. Cross 
section is 202,500 m 2 . 



layer 2 about 58 m 3 /d; it approaches steady state through 
all sections after 100 years with flux equaling about 
62.6 m 3 /d. The 100-year flux is 61.5 m 3 /d higher than 
the pre-injection flux because of the changed shale 
properties. 

Discussion 

The interpretative modeling completed herein has 
revealed several facts about fracking. First, MODFLOW 
can be coded to adequately simulate fracking. Simulated 
pressures are high, but velocities even near the well do 
not violate the assumptions for Darcian flow. Second, 
injection for 5 d causes extremely high pressure within 
the shale. The pressure decreases with distance from the 
well. The time to maximum pressure away from the well 
lags the time of maximum pressure at the well. The 
pressure drops back to close to its pre-injection level 
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at the well within 300 d, indicating the injection affects 
the flow for significantly longer periods than just during 
the fracking operation. Although the times may vary 
based on media properties, the difference would be at 
most a month or so, based on the various combinations 
of properties simulated. The system transitions within 
6 years due to changes in the shale properties. The 
equilibrium transport rate would transition from a system 
requiring thousands of years to one requiring less than 
100 years within less than 10 years. 

Third, most of the injected water in the simulation 
flows vertically rather than horizontally through the shale. 
This reflects the higher K ss 20 m above the well and the no 
flow boundary within 225 m laterally from the well, which 
emulates in situ shale properties that would manifest at 
some distance in the shale. 

Fourth, the interpretative model accurately and real- 
istically simulates long-term steady state flow conditions, 
with an upward flow that would advect whatever conser- 
vative constituents exist at depth. Using low, unfractured 
K values, the transport simulation may correspond with 
advective transport over geologic time although there are 
conditions for which it would occur much more quickly 
(Figure 4). If the is 0.01 m/d, transport could occur 
on the order of a few hundreds of years. Faults through the 
overburden could speed the transport time considerably. 
Reasonable scenarios presented herein suggest the travel 
time could be decreased further by an order of magnitude. 

Fifth, fracking increases the by several orders 
of magnitude. Out-of-formation fracking (Fisher and 
Warpinski 2011) would increase the K in the overburden, 
thereby changing the regional hydrogeology. Vertical flow 
could change over broad areas if the expected density 
of wells in the Marcellus shale region (NYDEC 2009) 
actually occurs. 

Sixth, if newly fractured shale or out-of-formation 
fractures come close to contacting fault fracture zones, 
contaminants could reach surface areas in tens of years, 
or less. Faults can decrease the simulated particle travel 
time several orders of magnitude. 



Conclusion 

Fracking can release fluids and contaminants from 
the shale either by changing the shale and overburden 
hydrogeology or simply by the injected fluid forcing other 
fluids out of the shale. The complexities of contaminant 
transport from hydraulically fractured shale to near- 
surface aquifers render estimates uncertain, but a range 
of interpretative simulations suggest that transport times 
could be decreased from geologic time scales to as 
few as tens of years. Preferential flow through natural 
fractures fracking-induced fractures could further decrease 
the travel times to as little as just a few years. 

There is no data to verify either the pre- or 
post-fracking properties of the shale. The evidence for 
potential vertical contaminant flow is strong, but there 
are also almost no monitoring systems that would 



detect contaminant transport as considered herein. Several 
improvements could be made. 

• Prior to hydraulic fracturing operations, the subsurface 
should be mapped for the presence of faults and 
measurement of their properties. 

• A reasonable setback distance from the fracking to 
the faults should be established. The setback distance 
should be based on a reasonable risk analysis of fracking 
increasing the pressures within the fault. 

• The properties of the shale should be verified, post- 
fracking, to assess how the hydrogeology will change. 

• A system of deep and shallow monitoring wells and 
piezometers should be established in areas expect- 
ing significant development, before that development 
begins (Williams 2010). 
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